#include"TShieldGenerator.h"

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// General Grid Tool Auxiliary Functions

void GridToolsGetStrides(int dim, int *dims, int *strides) {
	int d;
	strides[dim-1]=1;
	for(d=dim-2;d>=0;d--) {
		strides[d]=strides[d+1]*dims[d+1];
	}
}

int GridToolsGetIdFromPos(int dim, int *pos, int *strides) {
	int d, result;
	result=0;
	for(d=0;d<dim;d++) {
		result+=pos[d]*strides[d];
	}
	return result;
}

void GridToolsGetPosFromId(int dim, int id, int *pos, int *strides) {
	int d;
	pos[0]=id/strides[0];
	for(d=1;d<dim;d++) {
		pos[d]=(id%strides[d-1])/strides[d];
	}
}

//////////////////////////////

void GridToolsGetNeighbours(int dim, int *dims, TVarListHandler *neighbours) {
		GridToolsGetNeighbours_Torus(dim,dims,0,neighbours);
}

void GridToolsGetNeighbours_Torus(int dim, int *dims, int torusDim, TVarListHandler *neighbours) {
	/* dim: dimension of grid, dims: points per axis, torusDim: dimensions 0 to torusDim-1 are considered cyclical. */
	int *strides;
	strides=(int*) malloc(sizeof(int)*dim);
	int *pos;
	pos=(int*) malloc(sizeof(int)*dim);

	GridToolsGetStrides(dim,dims,strides);

	GridToolsGetNeighbours_Torus_iterateXVariables(neighbours, pos, dims, strides, dim, torusDim, 0);

	free(strides);
	free(pos);
}

void GridToolsGetNeighbours_Torus_iterateXVariables(TVarListHandler *neighbours, int *pos, int *dims, int *strides,
		int dim, int torusDim, int d) {
	int x;
	if(d<dim) {
		// as long as we have not reached the last axis, launch iteration along next axis
		for(x=0;x<dims[d];x++) {
			pos[d]=x;
			GridToolsGetNeighbours_Torus_iterateXVariables(neighbours,pos,dims,strides,dim,torusDim,d+1);
		}
	} else {
		// once a specific grid point is reached, add points
		//addVariables_Shields(xVars,xMap,xPos);
		int xId,d;
		xId=GridToolsGetIdFromPos(dim,pos,strides);
		for(d=0;d<dim;d++) {
			if(pos[d]>0) {
				neighbours->addToLine(xId,xId-strides[d]);
			} else {
				if(d<torusDim) {
					neighbours->addToLine(xId,xId+strides[d]*(dims[d]-1));
				}
			}
			if(pos[d]<dims[d]-1) {
				neighbours->addToLine(xId,xId+strides[d]);
			} else {
				if(d<torusDim) {
					neighbours->addToLine(xId,xId-strides[d]*(dims[d]-1));
				}
			}
		}

	}
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorBase */

void TShieldGeneratorBase::getXMap(int *xMap, TVarListHandler *xSupport) {
	int x;
	for(x=0;x<xSupport->res;x++) {
		xMap[x]=xSupport->varList[x]->at(0);
	}
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorGrid_SqrEuclidean */



TShieldGeneratorGrid_SqrEuclidean::TShieldGeneratorGrid_SqrEuclidean(int _dim, int *_xDims, int *_yDims) {
	dim=_dim;
	xDims=_xDims;
	yDims=_yDims;
	xStrides=(int*) malloc(sizeof(int)*dim);
	yStrides=(int*) malloc(sizeof(int)*dim);

	GridToolsGetStrides(dim,xDims,xStrides);
	GridToolsGetStrides(dim,yDims,yStrides);
}

TShieldGeneratorGrid_SqrEuclidean::~TShieldGeneratorGrid_SqrEuclidean() {
	free(xStrides);
	free(yStrides);
}

void TShieldGeneratorGrid_SqrEuclidean::generateShield(TVarListHandler *xVars, TVarListHandler *xSupport) {
	int *xPos,d;
	int *xMap;
	// A shield is generated by iterating over all x-variables.
	// For convenience this iteration is not simply carried out 1..xres
	// But in grid coordinates directly. Then the existence, location and number of neighbours can be
	// determined more directly.


	/* extract map from support */
	xMap=(int*) malloc(sizeof(int)*xSupport->res);
	getXMap(xMap,xSupport);

	// xPos will keep the current coordinates on the grid
	xPos=(int*) malloc(sizeof(int)*dim);
	for(d=0;d<dim;d++) {
		xPos[d]=0;
	}

	// start iterating along the first axis
	iterateXVariables(xVars, xMap, xPos, 0);

	free(xMap);
	free(xPos);

}

void TShieldGeneratorGrid_SqrEuclidean::iterateXVariables(TVarListHandler *xVars, int *xMap, int *xPos, int d) {
	int x;
	if(d<dim) {
		// as long as we have not reached the last axis, launch iteration along next axis
		for(x=0;x<xDims[d];x++) {
			xPos[d]=x;
			iterateXVariables(xVars,xMap,xPos,d+1);
		}
	} else {
		// once a specific grid point is reached, add points
		addVariables_Shields(xVars,xMap,xPos);
		addVariables_Rectangles(xVars,xMap,xPos);
	}
}

void TShieldGeneratorGrid_SqrEuclidean::addVariables_Shields(TVarListHandler *xVars, int *xMap, int *xPos) {
	int xId,d;
	xId=GridToolsGetIdFromPos(dim,xPos,xStrides);
	for(d=0;d<dim;d++) {
		if(xPos[d]>0) {
			xVars->addToLine(xId,xMap[xId-xStrides[d]]);
		}
		if(xPos[d]<xDims[d]-1) {
			xVars->addToLine(xId,xMap[xId+xStrides[d]]);
		}
	}
}

void TShieldGeneratorGrid_SqrEuclidean::addVariables_Rectangles(TVarListHandler *xVars, int *xMap, int *xPos) {
	int *yPos,*yMin,*yMax,d,yId,xId;

	yPos=(int*) malloc(sizeof(int)*dim);
	yMin=(int*) malloc(sizeof(int)*dim);
	yMax=(int*) malloc(sizeof(int)*dim);
	for(d=0;d<dim;d++) {
		yPos[d]=0;
		yMin[d]=0;
		yMax[d]=0;
	}

	// determine dimensions of cube on Y grid
	// first determine id of current xPos
	xId=GridToolsGetIdFromPos(dim,xPos,xStrides);
	// go through all dimensions
	for(d=0;d<dim;d++) {
		// along each dimension, check whether we're on lower boundary
		if(xPos[d]>0) {
			// if not, determine yId of lower neighbour
			yId=xMap[xId-xStrides[d]];
			// compute Pos variable
			GridToolsGetPosFromId(dim,yId,yPos,yStrides);
			// set yMin[d] to corresponding entry
			yMin[d]=yPos[d];
		} else {
			// if x is at lower boundary, all y need to be added
			yMin[d]=0;
		}
		// now same procedure for upper boundary
		if(xPos[d]<xDims[d]-1) {
			yId=xMap[xId+xStrides[d]];
			GridToolsGetPosFromId(dim,yId,yPos,yStrides);
			// add +1 here for later for loop compatibility
			yMax[d]=yPos[d]+1;
		} else {
			yMax[d]=yDims[d];
		}
	}

	// start iterating along the first axis
	iterateYVariables(xVars, xId,yPos,yMin,yMax, 0);

	free(yPos);
	free(yMin);
	free(yMax);
}


void TShieldGeneratorGrid_SqrEuclidean::iterateYVariables(TVarListHandler *xVars, int xId,
		int *yPos, int *yMin, int *yMax, int d) {
	int y;
	if(d<dim) {
		// as long as we have not reached the last axis, launch iteration along next axis
		for(y=yMin[d];y<yMax[d];y++) {
			yPos[d]=y;
			iterateYVariables(xVars,xId,yPos,yMin,yMax,d+1);
		}
	} else {
		int yId;
		yId=GridToolsGetIdFromPos(dim,yPos,yStrides);
		//cout << yId << " ";
		xVars->addToLine(xId,yId);

	}
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorGrid_Torus */

/*

TShieldGeneratorGrid_Torus::TShieldGeneratorGrid_Torus(int _dim, int *_xDims, int *_yDims, double *_radii) {
	dim=_dim;
	xDims=_xDims;
	yDims=_yDims;
	xStrides=(int*) malloc(sizeof(int)*dim);
	yStrides=(int*) malloc(sizeof(int)*dim);

	radii=_radii;

	GridToolsGetStrides(dim,xDims,xStrides);
	GridToolsGetStrides(dim,yDims,yStrides);
}

TShieldGeneratorGrid_Torus::~TShieldGeneratorGrid_Torus() {
	free(xStrides);
	free(yStrides);
}

void TShieldGeneratorGrid_Torus::generateShield(TVarListHandler *xVars, TVarListHandler *xSupport) {
	int *xPos,d;
	int *xMap;
	// A shield is generated by iterating over all x-variables.
	// For convenience this iteration is not simply carried out 1..xres
	// But in grid coordinates directly. Then the existence, location and number of neighbours can be
	// determined more directly.


	// extract map from support
	xMap=(int*) malloc(sizeof(int)*xSupport->res);
	getXMap(xMap,xSupport);

	// xPos will keep the current coordinates on the grid
	xPos=(int*) malloc(sizeof(int)*dim);
	for(d=0;d<dim;d++) {
		xPos[d]=0;
	}

	// start iterating along the first axis
	iterateXVariables(xVars, xMap, xPos, 0);

	free(xMap);
	free(xPos);

}

void TShieldGeneratorGrid_Torus::iterateXVariables(TVarListHandler *xVars, int *xMap, int *xPos, int d) {
	int x;
	if(d<dim) {
		// as long as we have not reached the last axis, launch iteration along next axis
		for(x=0;x<xDims[d];x++) {
			xPos[d]=x;
			iterateXVariables(xVars,xMap,xPos,d+1);
		}
	} else {
		// once a specific grid point is reached, add points
		addVariables_Shields(xVars,xMap,xPos);
		addVariables_Rectangles(xVars,xMap,xPos);
	}
}

void TShieldGeneratorGrid_Torus::addVariables_Shields(TVarListHandler *xVars, int *xMap, int *xPos) {
	int xId,d;
	xId=GridToolsGetIdFromPos(dim,xPos,xStrides);
	for(d=0;d<dim;d++) {
		if(xPos[d]>0) {
			xVars->addToLine(xId,xMap[xId-xStrides[d]]);
		}
		if(xPos[d]<xDims[d]-1) {
			xVars->addToLine(xId,xMap[xId+xStrides[d]]);
		}
	}
}

void TShieldGeneratorGrid_Torus::addVariables_Rectangles(TVarListHandler *xVars, int *xMap, int *xPos) {
	int *yPos,*yMin,*yMax,d,yId,xId;

	yPos=(int*) malloc(sizeof(int)*dim);
	yMin=(int*) malloc(sizeof(int)*dim);
	yMax=(int*) malloc(sizeof(int)*dim);
	for(d=0;d<dim;d++) {
		yPos[d]=0;
		yMin[d]=0;
		yMax[d]=0;
	}

	// determine dimensions of cube on Y grid
	// first determine id of current xPos
	xId=GridToolsGetIdFromPos(dim,xPos,xStrides);
	// go through all dimensions
	for(d=0;d<dim;d++) {
		// along each dimension, check whether we're on lower boundary
		if(xPos[d]>0) {
			// if not, determine yId of lower neighbour
			yId=xMap[xId-xStrides[d]];
			// compute Pos variable
			GridToolsGetPosFromId(dim,yId,yPos,yStrides);
			// set yMin[d] to corresponding entry
			yMin[d]=yPos[d];
		} else {
			// if x is at lower boundary, all y need to be added
			yMin[d]=0;
		}
		// now same procedure for upper boundary
		if(xPos[d]<xDims[d]-1) {
			yId=xMap[xId+xStrides[d]];
			GridToolsGetPosFromId(dim,yId,yPos,yStrides);
			// add +1 here for later for loop compatibility
			yMax[d]=yPos[d]+1;
		} else {
			yMax[d]=yDims[d];
		}
	}

	// start iterating along the first axis
	iterateYVariables(xVars, xId,yPos,yMin,yMax, 0);

	free(yPos);
	free(yMin);
	free(yMax);
}


void TShieldGeneratorGrid_Torus::iterateYVariables(TVarListHandler *xVars, int xId,
		int *yPos, int *yMin, int *yMax, int d) {
	int y;
	if(d<dim) {
		// as long as we have not reached the last axis, launch iteration along next axis
		for(y=yMin[d];y<yMax[d];y++) {
			yPos[d]=y;
			iterateYVariables(xVars,xId,yPos,yMin,yMax,d+1);
		}
	} else {
		int yId;
		yId=GridToolsGetIdFromPos(dim,yPos,yStrides);
		//cout << yId << " ";
		xVars->addToLine(xId,yId);

	}
}
*/

/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorGrid_Padding */



TShieldGeneratorGrid_Padding::TShieldGeneratorGrid_Padding(int _dim, int *_xDims, int *_yDims, int _width) {
	dim=_dim;
	xDims=_xDims;
	yDims=_yDims;
	width=_width;
	xStrides=(int*) malloc(sizeof(int)*dim);
	yStrides=(int*) malloc(sizeof(int)*dim);

	GridToolsGetStrides(dim,xDims,xStrides);
	GridToolsGetStrides(dim,yDims,yStrides);
}

TShieldGeneratorGrid_Padding::~TShieldGeneratorGrid_Padding() {
	free(xStrides);
	free(yStrides);
}

void TShieldGeneratorGrid_Padding::generateShield(TVarListHandler *xVars,
		__attribute__((unused)) TVarListHandler *xSupport) {

	TVarListHandler *xSupportOld; // store copy of initial support
	int n; // count loops until width

	n=0;
	while(n<width) {
		// copy support to avoid self-interference of iteration while adding variables
		xSupportOld=new TVarListHandler(xVars);

		// iterate over variables
		iterateVariables(xVars, xSupportOld);

		delete(xSupportOld);
		n++;
	}
}

void TShieldGeneratorGrid_Padding::iterateVariables(TVarListHandler *xVars,
		TVarListHandler *xSupport) {

	int *xPos,*yPos;
	int x,y,yIndex;

	xPos=(int*) malloc(sizeof(int)*dim);
	yPos=(int*) malloc(sizeof(int)*dim);

	for(x=0;x < (int) xVars->res;x++) {
		// iterate over x variables

		// determine xPos
		GridToolsGetPosFromId(dim, x, xPos, xStrides);

		for(yIndex=0;yIndex < (int) xSupport->lenList->at(x);yIndex++) {
			// iterate over row
			y=xSupport->varList[x]->at(yIndex);

			// determine yPos
			GridToolsGetPosFromId(dim, y, yPos, yStrides);

			addVariables(xVars,x,y,xPos,yPos);
		}
	}

	free(xPos);
	free(yPos);
}

void TShieldGeneratorGrid_Padding::addVariables(TVarListHandler *xVars,
		int xId, int yId,
		int *xPos, int *yPos) {

	int d;


	for(d=0;d<dim;d++) {

		// for support element described by xPos, yPos, add x-neighbour elements
		if(xPos[d]>0) {
			xVars->addToLine(xId-xStrides[d],yId);
		}
		if(xPos[d]<xDims[d]-1) {
			xVars->addToLine(xId+xStrides[d],yId);
		}


		// for support element described by xPos, yPos, add y-neighbour elements
		if(yPos[d]>0) {
			xVars->addToLine(xId,yId-yStrides[d]);
		}
		if(yPos[d]<yDims[d]-1) {
			xVars->addToLine(xId,yId+yStrides[d]);
		}

		// neighbour-neighbour elements
		/*
		int d2;
		if(xPos[d]>0) {
			for(d2=0;d2<dim;d2++) {
				if(yPos[d2]>0) {
					xVars->addToLine(xId-xStrides[d],yId-yStrides[d2]);
				}
				if(yPos[d2]<yDims[d2]-1) {
					xVars->addToLine(xId-xStrides[d],yId+yStrides[d2]);
				}
			}
		}
		if(xPos[d]<xDims[d]-1) {
			for(d2=0;d2<dim;d2++) {
				if(yPos[d2]>0) {
					xVars->addToLine(xId+xStrides[d],yId-yStrides[d2]);
				}
				if(yPos[d2]<yDims[d2]-1) {
					xVars->addToLine(xId+xStrides[d],yId+yStrides[d2]);
				}
			}
		}
		*/

	}
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorTreeBase */


TShieldGeneratorTreeBase::TShieldGeneratorTreeBase(int _dim,
		THierarchicalPartition* _yPartition, double** _yPos, double** _yRadii,
		int _lBottom, int _lTop, double* _xpos, TVarListHandler* _xNeighbours) {

	dim=_dim;
	yPartition=_yPartition;
	yPos=_yPos;
	yRadii=_yRadii;
	lBottom=_lBottom;
	lTop=_lTop;
	xpos=_xpos;
	xNeighbours=_xNeighbours;
	xres=xNeighbours->res;

}

TShieldGeneratorTreeBase::~TShieldGeneratorTreeBase() {
}

void TShieldGeneratorTreeBase::generateShield(TVarListHandler *xVars, TVarListHandler *xSupport)
 {
	int x;
	int *xMap;

	/* extract map from support */
	xMap=(int*) malloc(sizeof(int)*xSupport->res);
	getXMap(xMap,xSupport);


	for(x=0;x<xres;x++) {
		addVariables_Shields(xVars, xMap, x);
		addVariables_Polytopes(xVars, xMap, x);
	}

	free(xMap);

}

void TShieldGeneratorTreeBase::addVariables_Shields(
		TVarListHandler* xVars, int* xMap, int x) {
	int nIndex,xn;
	// iterate over neighbours of x
	for(nIndex=0;nIndex<xNeighbours->lenList->at(x);nIndex++) {
		// extract id of neighbour
		xn=xNeighbours->varList[x]->at(nIndex);
		//cout << "\t\t" << xn << endl;
		xVars->addToLine(x,xMap[xn]);
	}
}

void TShieldGeneratorTreeBase::addVariables_Polytopes(
		TVarListHandler* xVars, int* xMap, int x) {
	int yBParent;
	for(yBParent=0;yBParent<yPartition->layers[lTop]->nCells;yBParent++) {
		iterateYVariables(xVars,xMap,x,lTop,yBParent);
	}
}

void TShieldGeneratorTreeBase::iterateYVariables(TVarListHandler* xVars,
		int* xMap, int xA, int lParent, int yBParent) {

	int yB,yBIndex;
	// this function is to be called only for coarser layers, not on finest layer.

	// for convenience, introduce direct references to relevant layer
	TPartitionLayer *parentLayer;
	parentLayer=yPartition->layers[lParent];
	// go through all children of current node
	for(yBIndex=0;yBIndex<parentLayer->nChildren[yBParent];yBIndex++) {
		// extract id
		yB=parentLayer->children[yBParent][yBIndex];
		// call hierarchical shielding condition test
		if(!checkCondition(xA,lParent+1,yB,xMap)) {
			// if condition fails
			if(lParent<lBottom-1) {
				// do test at finer scale
				iterateYVariables(xVars,xMap,xA,lParent+1,yB);
			} else {
				// when on finest scale, add variable
				xVars->addToLine(xA,yB);
			}
		}
	}

}


bool TShieldGeneratorTreeBase::checkCondition(int xA, int l, int yB,
		int* xMap) {

	int nIndex,xn,yn;

	//nCalls_checkCondition[l-lTop]++;

	for(nIndex=0;nIndex<xNeighbours->lenList->at(xA);nIndex++) {
		// extract id of neighbour
		xn=xNeighbours->varList[xA]->at(nIndex);
		yn=xMap[xn];
		//nCalls_checkConditionPlane[l-lTop]++;
		if(checkConditionPlane(xA,xn,l,yB,yn)) {
			//cout << "\t\t(" << l << "," << yB << ") : true" << endl;
			return true;
		}
	}

	//cout << "\t\t(" << l << "," << yB << ") : false" << endl;
	return false;
}

bool TShieldGeneratorTreeBase::checkConditionPlane(__attribute__((unused)) int xA,
		__attribute__((unused)) int x, __attribute__((unused)) int l,
		__attribute__((unused)) int yB, __attribute__((unused)) int y) {
	return false;
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorTree_SqrEuclidean */

template<class TReferenceTree>
TShieldGeneratorTree_SqrEuclideanPrototype<TReferenceTree>::TShieldGeneratorTree_SqrEuclideanPrototype(int _dim,
		THierarchicalPartition* _yPartition, double** _yPos, double** _yRadii,
		int _lBottom, int _lTop, double* _xpos, TVarListHandler* _xNeighbours) :
			TReferenceTree(_dim, _yPartition, _yPos, _yRadii,
						_lBottom, _lTop, _xpos, _xNeighbours) {
}

template<class TReferenceTree>
bool TShieldGeneratorTree_SqrEuclideanPrototype<TReferenceTree>::checkConditionPlane(int xA, int x, int l,
		int yB, int y) {
	// xA: index of initial x
	// yB: index of final y
	// x,y: indices of shielding pair
	// l: hierarchy level. the higher, the finer the resolution.
	double dist,norm,xpart;
	int i;
	// xpart: x-xA (called a in notes)
	// dist: (yB-y).(x-xA)
	// norm: in the end: length of x-xA (sqrt only taken in final else statement)
	dist=0;
	norm=0;
	for(i=0;i<dim;i++) {
		xpart=(xpos[x*dim+i]-xpos[xA*dim+i]);
		dist+=xpart*(yPos[l][yB*dim+i]-yPos[lBottom][y*dim+i]);
		norm+=xpart*xpart;
	}
	if(l>=lBottom) {
		return (dist>shieldingTolerance);
	} else {
		norm=sqrt(norm);
		return ((dist-norm*yRadii[l][yB])>shieldingTolerance);
	}
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorTree_Torus */

template<class TReferenceTree>
TShieldGeneratorTree_TorusPrototype<TReferenceTree>::TShieldGeneratorTree_TorusPrototype(int _dim,
		THierarchicalPartition* _yPartition, double** _yPos, double** _yRadii,
		double *** _yTorusRadii,
		int _lBottom, int _lTop, double* _xpos, TVarListHandler* _xNeighbours,
		double *_torusRadii, int _torusDim) :
			TReferenceTree(_dim, _yPartition, _yPos, _yRadii,
						_lBottom, _lTop, _xpos, _xNeighbours) {
	yTorusRadii=_yTorusRadii;
	torusRadii=_torusRadii;
	torusDim=_torusDim;
}

template<class TReferenceTree>
bool TShieldGeneratorTree_TorusPrototype<TReferenceTree>::checkCondition(int xA, int l, int yB, int *xMap) {
	int axis,nIndex,xn,yn;
	double dist;

	//cout << "\t(" << xA << "," << l << "," << yB << ")" << endl;

	for(nIndex=0;nIndex<xNeighbours->lenList->at(xA);nIndex++) {
		// reset slack
		dist=0;

		// extract id of neighbour
		xn=xNeighbours->varList[xA]->at(nIndex);
		yn=xMap[xn];

		//cout << "\t\t(" << xA << "," << xn << "," << l << "," << yB << "," << yn << ")" << endl;
		// get psi slack of torus-y-fied dimensions.
		for(axis=0;axis<torusDim;axis++) {
			dist+=slackConditionS1(xA,xn,l,yB,yn,axis);
		}
		if(torusDim<dim) {
			dist+=slackConditionPlane(xA,xn,l,yB,yn);
		}
		if(dist>shieldingTolerance) {
			return true;
		}
	}
	return false;
}

template<class TReferenceTree>
double TShieldGeneratorTree_TorusPrototype<TReferenceTree>::slackConditionPlane(int xA, int x, int l,
		int yB, int y) {
	// basically a copy of SqrEuclidean::checkConditionPlane, just works on the last dimensions, that are not
	//   torus-y-fied
	// xA: index of initial x
	// yB: index of final y
	// x,y: indices of shielding pair
	// l: hierarchy level. the higher, the finer the resolution.
	double dist,norm,xpart;
	int i;
	// xpart: x-xA (called a in notes)
	// dist: (yB-y).(x-xA)
	// norm: in the end: length of x-xA (sqrt only taken in final else statement)
	dist=0;
	norm=0;
	for(i=torusDim;i<dim;i++) {
		xpart=(xpos[x*dim+i]-xpos[xA*dim+i]);
		dist+=xpart*(yPos[l][yB*dim+i]-yPos[lBottom][y*dim+i]);
		norm+=xpart*xpart;
	}
//	if(l>=lBottom) {
//		return (dist>shieldingTolerance);
//	} else {
//		norm=sqrt(norm);
//		return ((dist-norm*yRadii[l][yB])>shieldingTolerance);
//	}
	if(l<lBottom) {
		norm=sqrt(norm);
		dist=dist-norm*yRadii[l][yB];
	}
	//cout << "\t\t\tEuclidean: " << dist << endl;
	return dist;
}


template<class TReferenceTree>
double TShieldGeneratorTree_TorusPrototype<TReferenceTree>::slackConditionS1(int xA, int x, int l,
		int yB, int y, int axis) {
	// do a S1-type shielding check with (xA,x,yB,y). checkCondition will call this function such that
	//   x is a neighbour along axis, so that the S1 formulas can be applied.
	// xA: index of initial x
	// yB: index of final y
	// x,y: indices of shielding pair
	// l: hierarchy level. the higher, the finer the resolution.
	double delta,coord_xA,coord_x,coord_y,coord_yB,dist;

	coord_xA=xpos[xA*dim+axis];
	coord_x=xpos[x*dim+axis];
	coord_y=yPos[lBottom][y*dim+axis];
	coord_yB=yPos[l][yB*dim+axis];

	//cout << "\t\t\t" << "(" << coord_xA << "," << coord_x << "," << coord_y << "," << coord_yB << ")" << endl;

	// transformation into standard coordinate system
	// range=[0,1], xA=0.5, x>=0.5
	delta=0.5-coord_xA/torusRadii[axis];
	coord_xA=0.5;
	coord_x=translateCoord(coord_x,torusRadii[axis],delta);
	coord_y=translateCoord(coord_y,torusRadii[axis],delta);
	coord_yB=translateCoord(coord_yB,torusRadii[axis],delta);
	if(coord_x<0.5) {
		coord_x=1-coord_x;
		coord_y=1-coord_y;
		coord_yB=1-coord_yB;
	}

	//cout << "\t\t\t" << "(" << coord_xA << "," << coord_x << "," << coord_y << "," << coord_yB << ")" << endl;

	double psi_y,psi_yB;
	psi_y=psi(coord_y,coord_x,0);

	if(l>=lBottom) {
		psi_yB=psi(coord_yB,coord_x,0);
	} else {
		psi_yB=psi(coord_yB,coord_x,yTorusRadii[axis][l][yB]/torusRadii[axis]);
	}

	//cout << "\t\t\t" << "(" << psi_y << "," << psi_yB << ")" << endl;
	dist=psi_yB-psi_y;
	//return (dist>shieldingTolerance);

	// must not forget to rescale dist, to make it comparable with other axes.
	return dist*pow(torusRadii[axis],2);
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorTree_SqrEuclideanNoise */


template<class TReferenceTree>
TShieldGeneratorTree_SqrEuclideanNoisePrototype<TReferenceTree>::TShieldGeneratorTree_SqrEuclideanNoisePrototype(int _dim,
		THierarchicalPartition* _yPartition, double** _yPos, double** _yRadii,
		int _lBottom, int _lTop, double* _xpos, TVarListHandler* _xNeighbours,
		double **_c, double _eta, double _lambda) :
			TReferenceTree(_dim, _yPartition, _yPos, _yRadii,
						_lBottom, _lTop, _xpos, _xNeighbours) {
	c=_c;
	eta=_eta;
	lambda=_lambda;
	yres=_yPartition->layers[_lBottom]->nCells;
}

template<class TReferenceTree>
bool TShieldGeneratorTree_SqrEuclideanNoisePrototype<TReferenceTree>::checkConditionPlane(int xA, int x, int l,
		int yB, int y) {
	// xA: index of initial x
	// yB: index of final y
	// x,y: indices of shielding pair
	// l: hierarchy level. the higher, the finer the resolution.
	double dist,norm,xpart;
	int i;
	// xpart: x-xA (called a in notes)
	// dist: (yB-y).(x-xA)
	// norm: in the end: length of x-xA (sqrt only taken in final else statement)
	//EUCL_lincomb(xpos+(x*dim),xpos+(xA*dim),xpart,1.,-1.,dim);


	if(l>=lBottom) {
		dist=c[l][xA*yres+yB]-c[l][x*yres+yB]-c[l][xA*yres+y]+c[l][x*yres+y];
		return (dist>shieldingTolerance);
	} else {
		dist=0;
		norm=0;
		for(i=0;i<dim;i++) {
			xpart=(xpos[x*dim+i]-xpos[xA*dim+i]);
			dist+=xpart*(yPos[l][yB*dim+i]-yPos[lBottom][y*dim+i]);
			norm+=xpart*xpart;
		}

		norm=sqrt(norm);
		return ((dist-norm*yRadii[l][yB])-eta-lambda*norm>shieldingTolerance);
	}
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorTree_PEuclidean */

template<class TReferenceTree>
TShieldGeneratorTree_PEuclideanPrototype<TReferenceTree>::TShieldGeneratorTree_PEuclideanPrototype(int _dim,
		THierarchicalPartition* _yPartition, double** _yPos, double** _yRadii,
		int _lBottom, int _lTop, double* _xpos, TVarListHandler* _xNeighbours,
		double _p, double _slack) :
				TReferenceTree(_dim, _yPartition, _yPos, _yRadii,
						_lBottom, _lTop, _xpos, _xNeighbours) {

	p=_p;
	slack=_slack;
	xAyB=(double*) malloc(sizeof(double)*dim);
	xAys=(double*) malloc(sizeof(double)*dim);
	xsyB=(double*) malloc(sizeof(double)*dim);
	xsys=(double*) malloc(sizeof(double)*dim);
	xAxs=(double*) malloc(sizeof(double)*dim);
}

template<class TReferenceTree>
TShieldGeneratorTree_PEuclideanPrototype<TReferenceTree>::~TShieldGeneratorTree_PEuclideanPrototype() {
	free(xAyB);
	free(xAys);
	free(xsyB);
	free(xsys);
	free(xAxs);
}

template<class TReferenceTree>
bool TShieldGeneratorTree_PEuclideanPrototype<TReferenceTree>::checkConditionPlane(int xA, int x, int l,
		int yB, int y) {
	// xA: index of initial x
	// yB: index of final y
	// x,y: indices of shielding pair
	// l: hierarchy level. the higher, the finer the resolution.
	double dist;


	// compute necessary vector differences
	TReferenceTree::EUCL_lincomb(xpos+(xA*dim),yPos[lBottom]+(y*dim),xAys,1.,-1.,dim); // xA-ys
	TReferenceTree::EUCL_lincomb(xpos+(x*dim),yPos[lBottom]+(y*dim),xsys,1.,-1.,dim); // xs-ys
	TReferenceTree::EUCL_lincomb(xpos+(x*dim),yPos[l]+(yB*dim),xsyB,1.,-1.,dim); // xs-yB

	// compute exact r.h.s. of shielding condition
	dist=-(getPhi(xAys)-getPhi(xsys));

	if(l>=lBottom) {
		// if at finest level,

		// compute necessary vector differences (additional)
		TReferenceTree::EUCL_lincomb(xpos+(xA*dim),yPos[l]+(yB*dim),xAyB,1.,-1.,dim); // xA-yB

		// exact l.h.s. of shielding condition
		dist+=(getPhi(xAyB)-getPhi(xsyB));
		return (dist>=shieldingTolerance+slack);

	} else {

		// if at higher levels

		// compute necessary vector differences (additional)
		TReferenceTree::EUCL_lincomb(xpos+(xA*dim),xpos+(x*dim),xAxs,1.,-1.,dim); // xA-xs


		double xi,A,B,xsyB_Norm,xAxs_Norm,xAxs_xsyB_Angle;
		xsyB_Norm=sqrt(TReferenceTree::EUCL_innerProduct(xsyB,xsyB,dim));
		xAxs_Norm=sqrt(TReferenceTree::EUCL_innerProduct(xAxs,xAxs,dim));

		// compute xi (maximal angle distortion of xsyB by delta)
		if(yRadii[l][yB]<xsyB_Norm) {
			xi=asin(yRadii[l][yB]/xsyB_Norm);
		} else {
			xi=M_PI;
		}

		// compute angle between a and vB (set to 0, if one vector is zero)
		if((xAxs_Norm>0) && (xsyB_Norm>0)) {
			xAxs_xsyB_Angle=acos(TReferenceTree::EUCL_innerProduct(xAxs,xsyB,dim)/(xAxs_Norm*xsyB_Norm));
		} else {
			xAxs_xsyB_Angle=0.;
		}

		// compute factor A
		if(xAxs_xsyB_Angle+xi<M_PI) {
			A=xAxs_Norm*cos(xAxs_xsyB_Angle+xi);
		} else {
			A=-xAxs_Norm;
		}

		// compute factor B
		if(A>=0) {
			B=max(0.,xsyB_Norm-yRadii[l][yB]);
		} else {
			B=xsyB_Norm+yRadii[l][yB];
		}

		dist+=p*pow(B,p-1)*A;

		return (dist>shieldingTolerance+slack);

	}
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorTree_Sphere */


template<class TReferenceTree>
TShieldGeneratorTree_SpherePrototype<TReferenceTree>::TShieldGeneratorTree_SpherePrototype(int _dim,
		THierarchicalPartition* _yPartition, double** _yPos, double** _yRadii,
		int _lBottom, int _lTop, double* _xpos, TVarListHandler* _xNeighbours,
		double _p) :
				TReferenceTree(_dim, _yPartition, _yPos, _yRadii,
						_lBottom, _lTop, _xpos, _xNeighbours) {

	p=_p;
}

template<class TReferenceTree>
TShieldGeneratorTree_SpherePrototype<TReferenceTree>::~TShieldGeneratorTree_SpherePrototype() {
}


template<class TReferenceTree>
bool TShieldGeneratorTree_SpherePrototype<TReferenceTree>::checkConditionPlane(int xA, int x, int l,
		int yB, int y) {
	// xA: index of initial x
	// yB: index of final y
	// x,y: indices of shielding pair
	// l: hierarchy level. the higher, the finer the resolution.

	// inner products
	double IPxAy, IPxy, IPxAx, IPxAyB, IPxyB;
	// final shielding bound estimate
	double dist;

	IPxAy=TReferenceTree::EUCL_innerProduct(xpos+(xA*dim),yPos[lBottom]+(y*dim),dim);
	IPxy=TReferenceTree::EUCL_innerProduct(xpos+(x*dim),yPos[lBottom]+(y*dim),dim);


	IPxAx=TReferenceTree::EUCL_innerProduct(xpos+(xA*dim),xpos+(x*dim),dim);
	IPxAyB=TReferenceTree::EUCL_innerProduct(xpos+(xA*dim),yPos[l]+(yB*dim),dim);
	IPxyB=TReferenceTree::EUCL_innerProduct(xpos+(x*dim),yPos[l]+(yB*dim),dim);


	// compute r.h.s. of shielding condition
	dist=-(getPhi(acos(IPxAy))-getPhi(acos(IPxy)));

	if(l>=lBottom) {
		// if at finest level,
		dist+=getPhi(acos(IPxAyB))-getPhi(acos(IPxyB));
		return (dist>=shieldingTolerance);
	} else {

		double theta, phi;
		//double thetaS;
		double deltaTheta, deltaPhi;
		double subgrad;
		double thetaMin, phiMax;
		double deltaD;

		double cosDeltaTheta,sinThetaS,sinTheta;

		theta=acos(IPxAyB);
		//thetaS=acos(IPxAx);
		sinTheta=std::sqrt(1-IPxAyB*IPxAyB);
		sinThetaS=std::sqrt(1-IPxAx*IPxAx);

		deltaTheta=yRadii[l][yB];
		cosDeltaTheta=cos(deltaTheta);

		phi=acos(std::min(1.,std::max(-1.,(IPxyB-IPxAx*IPxAyB)/(sinTheta*sinThetaS ) ) ) );

		// compute delta phi


		double cosDeltaThetaSqr,cosThetaSqr;
		// if point lies close to one of the poles
		cosDeltaThetaSqr=cosDeltaTheta*cosDeltaTheta;
		cosThetaSqr=IPxAyB*IPxAyB;
		if(cosThetaSqr>=cosDeltaThetaSqr) {
			deltaPhi=M_PI;
		} else {
			deltaPhi=acos(sqrt((cosDeltaThetaSqr-cosThetaSqr)/(1-cosThetaSqr)));

		}


		thetaMin=std::max(0.,theta-deltaTheta);
		phiMax=std::min(M_PI,phi+deltaPhi);


		double d1,d2;
		d1=thetaMin;
		d2=acos(sinThetaS*sin(thetaMin)*cos(phiMax)+IPxAx*cos(thetaMin));

		/*
		// fixed square variant
		deltaD=d1*d1-d2*d2;
		dist+=deltaD;
		return (dist>shieldingTolerance);
		*/


		// general subgradient variant

		deltaD=(d1-d2);

		// get subgradient bound
		if(deltaD>0) {
			subgrad=getSubgrad(std::max(0.,acos(IPxyB)-deltaTheta));
		} else {
			subgrad=getSubgrad(std::min(M_PI,acos(IPxyB)+deltaTheta));
		}

		dist+=subgrad*deltaD;

		return (dist>shieldingTolerance);

	}
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* TShieldGeneratorTree_Reflector */


template<class TReferenceTree>
TShieldGeneratorTree_ReflectorPrototype<TReferenceTree>::TShieldGeneratorTree_ReflectorPrototype(int _dim,
		THierarchicalPartition* _yPartition, double** _yPos, double** _yRadii,
		int _lBottom, int _lTop, double* _xpos, TVarListHandler* _xNeighbours
		) :
				TReferenceTree(_dim, _yPartition, _yPos, _yRadii,
						_lBottom, _lTop, _xpos, _xNeighbours) {

}

template<class TReferenceTree>
TShieldGeneratorTree_ReflectorPrototype<TReferenceTree>::~TShieldGeneratorTree_ReflectorPrototype() {
}


template<class TReferenceTree>
bool TShieldGeneratorTree_ReflectorPrototype<TReferenceTree>::checkConditionPlane(int xA, int x, int l,
		int yB, int y) {
	// xA: index of initial x
	// yB: index of final y
	// x,y: indices of shielding pair
	// l: hierarchy level. the higher, the finer the resolution.

	// inner products
	double IPxAy, IPxy, IPxAx, IPxAyB, IPxyB;
	// final shielding bound estimate
	double dist;

	IPxAy=TReferenceTree::EUCL_innerProduct(xpos+(xA*dim),yPos[lBottom]+(y*dim),dim);
	IPxy=TReferenceTree::EUCL_innerProduct(xpos+(x*dim),yPos[lBottom]+(y*dim),dim);


	IPxAx=TReferenceTree::EUCL_innerProduct(xpos+(xA*dim),xpos+(x*dim),dim);
	IPxAyB=TReferenceTree::EUCL_innerProduct(xpos+(xA*dim),yPos[l]+(yB*dim),dim);
	IPxyB=TReferenceTree::EUCL_innerProduct(xpos+(x*dim),yPos[l]+(yB*dim),dim);


	// compute r.h.s. of shielding condition
	dist=-(getPhi(acos(IPxAy))-getPhi(acos(IPxy)));

	if(l>=lBottom) {
		// if at finest level,
		dist+=getPhi(acos(IPxAyB))-getPhi(acos(IPxyB));
		return (dist>=shieldingTolerance);
	} else {

		double theta, phi;
		//double thetaS;
		double deltaTheta, deltaPhi;
		double subgrad;
		double thetaMax, phiMin;
		double deltaD;

		double cosDeltaTheta,sinThetaS,sinTheta;

		theta=acos(IPxAyB);
		//thetaS=acos(IPxAx);
		sinTheta=std::sqrt(1-IPxAyB*IPxAyB);
		sinThetaS=std::sqrt(1-IPxAx*IPxAx);


		deltaTheta=yRadii[l][yB];
		cosDeltaTheta=cos(deltaTheta);

		phi=acos(std::min(1.,std::max(-1.,(IPxyB-IPxAx*IPxAyB)/(sinTheta*sinThetaS ) ) ) );

		// compute delte phi

		double cosDeltaThetaSqr,cosThetaSqr;
		// if point lies close to one of the poles
		cosDeltaThetaSqr=cosDeltaTheta*cosDeltaTheta;
		cosThetaSqr=IPxAyB*IPxAyB;
		if(cosThetaSqr>=cosDeltaThetaSqr) {
			deltaPhi=M_PI;
		} else {
			deltaPhi=acos(sqrt((cosDeltaThetaSqr-cosThetaSqr)/(1-cosThetaSqr)));

		}


		thetaMax=std::min(M_PI,theta+deltaTheta);
		phiMin=std::max(0.,phi-deltaPhi);


		double d1,d2;
		d1=thetaMax;
		d2=acos(sinThetaS*sin(thetaMax)*cos(phiMin)+IPxAx*cos(thetaMax));

		// general subgradient variant

		deltaD=(d1-d2);

		// get subgradient bound
		if(deltaD>0) {
			subgrad=getSubgrad(std::max(0.,acos(IPxyB)-deltaTheta));
		} else {
			subgrad=getSubgrad(std::min(M_PI,acos(IPxyB)+deltaTheta));
		}

		dist+=subgrad*deltaD;

		return (dist>shieldingTolerance);

	}
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* Auxiliary Shielding Generators */


/* Comparision */


TShieldGeneratorComparison::TShieldGeneratorComparison(TShieldGeneratorBase *_shieldGenerator1, TShieldGeneratorBase *_shieldGenerator2) {
		shieldGenerator1=_shieldGenerator1;
		shieldGenerator2=_shieldGenerator2;
}

void TShieldGeneratorComparison::generateShield(TVarListHandler *xVars, TVarListHandler *xSupport) {
	TVarListHandler *xVarsCopy;
	xVarsCopy=new TVarListHandler(xVars);
	TVarListHandler *xSupportCopy;
	xSupportCopy=new TVarListHandler(xSupport);
	shieldGenerator1->generateShield(xVars,xSupport);
	shieldGenerator2->generateShield(xVarsCopy,xSupportCopy);
	xVars->sort();
	xVarsCopy->sort();

	bool difference=false;
	int x,yIndex,len1,len2;
	for(x=0;x<xVars->res;x++) {
		len1=xVars->lenList->at(x);
		len2=xVarsCopy->lenList->at(x);
		if(len1!=len2) {
			cout << "x=" << x << " : lengths differ" << endl;
			difference=true;
		} else {
			for(yIndex=0;yIndex<len1;yIndex++) {
				if(xVars->varList[x]->at(yIndex)!=xVarsCopy->varList[x]->at(yIndex)) {
					cout << "x=" << x << " yIndex=" << yIndex << " : entries differ" << endl;
					difference=true;
				}
			}
		}
	}
	if(!difference) {
		cout << "no difference detected in shields" << endl;
	} else {
		cout << "difference detected" << endl;
	}
	delete xVarsCopy;
	delete xSupportCopy;
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* Verification */

TShieldingVerification::TShieldingVerification(double* _c, int _xres, int _yres,
		TVarListHandler* _xNeighbours) {
	c=_c;
	xres=_xres;
	yres=_yres;
	xNeighbours=_xNeighbours;
}

TVarListHandler* TShieldingVerification::verify(TVarListHandler* xVars, int *xMap) {

	// THIS METHOD ASSUMES THAT xVars IS SORTED!

	int xA, yB, xsIndex, xs, ys;
	bool foundShield;
	TVarListHandler *result;
	result=new TVarListHandler();
	result->setupEmpty(xres);

	// variables to perform search of yB in xVars[xA]
	int varPos,currentY;

	for(xA=0;xA<xres;xA++) {
		varPos=0;
		if((int) xVars->lenList->at(xA)>varPos) {
			currentY=xVars->varList[xA]->at(varPos);
		} else {
			currentY=-1;
		}
		for(yB=0;yB<yres;yB++) {
			if(currentY==yB) {
				// if yB is at current position in varList
				varPos++;
				if((int) xVars->lenList->at(xA)>varPos) {
					currentY=xVars->varList[xA]->at(varPos);
				} else {
					currentY=-1;
				}
			} else {
				// else, do shielding check
				foundShield=false;
				for(xsIndex=0;xsIndex<(int) xNeighbours->lenList->at(xA);xsIndex++) {
					xs=xNeighbours->varList[xA]->at(xsIndex);
					ys=xMap[xs];
					if(c[xA*yres+yB]-c[xs*yres+yB]>c[xA*yres+ys]-c[xs*yres+ys]) {
						foundShield=true;
					}
				}
				if(!foundShield) {
					result->addToLine(xA,yB,false);
				}
			}
		}
	}
	return result;
}


TShieldingVerificationDuplex::TShieldingVerificationDuplex(double *_c, int _xres, int _yres,
		TVarListHandler *_xNeighbours, TVarListHandler *_yNeighbours) {
	c=_c;
	xres=_xres;
	yres=_yres;
	xNeighbours=_xNeighbours;
	yNeighbours=_yNeighbours;
}

TVarListHandler* TShieldingVerificationDuplex::verify(TVarListHandler *xVars,
		TVarListHandler *xSupport, TVarListHandler *ySupport) {

	// THIS METHOD ASSUMES THAT xVars IS SORTED!

	int xA, yB, xsIndex, xs, ysIndex, ys;
	bool foundShield;
	TVarListHandler *result;
	result=new TVarListHandler();
	result->setupEmpty(xres);

	// variables to perform search of yB in xVars[xA]
	int varPos,currentY;

	for(xA=0;xA<xres;xA++) {
		varPos=0;
		if((int) xVars->lenList->at(xA)>varPos) {
			currentY=xVars->varList[xA]->at(varPos);
		} else {
			currentY=-1;
		}
		for(yB=0;yB<yres;yB++) {
			if(currentY==yB) {
				// if yB is at current position in varList
				varPos++;
				if((int) xVars->lenList->at(xA)>varPos) {
					currentY=xVars->varList[xA]->at(varPos);
				} else {
					currentY=-1;
				}
			} else {
				// else, do shielding check
				foundShield=false;

				// do "box in Y"-check
				for(xsIndex=0;(xsIndex<(int) xNeighbours->lenList->at(xA)) && (!foundShield);xsIndex++) {
					xs=xNeighbours->varList[xA]->at(xsIndex);
					for(ysIndex=0;(ysIndex<(int) xSupport->lenList->at(xs)) && (!foundShield);ysIndex++) {
						ys=xSupport->varList[xs]->at(ysIndex);
						if(c[xA*yres+yB]-c[xs*yres+yB]>c[xA*yres+ys]-c[xs*yres+ys]) {
							foundShield=true;
						}
					}
				}

				// do "box in X"-check
				for(ysIndex=0;(ysIndex<(int) yNeighbours->lenList->at(yB)) && (!foundShield);ysIndex++) {
					ys=yNeighbours->varList[yB]->at(ysIndex);
					for(xsIndex=0;(xsIndex<(int) ySupport->lenList->at(ys)) && (!foundShield);xsIndex++) {
						xs=ySupport->varList[ys]->at(xsIndex);
						if(c[xA*yres+yB]-c[xs*yres+yB]>c[xA*yres+ys]-c[xs*yres+ys]) {
							foundShield=true;
						}
					}
				}


				if(!foundShield) {
					result->addToLine(xA,yB,false);
				}
			}
		}
	}
	return result;
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* Tree Benchmark Base Class */


TShieldGeneratorTreeBase_Benchmark::TShieldGeneratorTreeBase_Benchmark(int _dim,
		THierarchicalPartition* _yPartition, double** _yPos, double** _yRadii,
		int _lBottom, int _lTop, double* _xpos, TVarListHandler* _xNeighbours) :
		TShieldGeneratorTreeBase(_dim,
				_yPartition, _yPos, _yRadii,
				_lBottom, _lTop, _xpos, _xNeighbours) {
}

TShieldGeneratorTreeBase_Benchmark::~TShieldGeneratorTreeBase_Benchmark() {
}

void TShieldGeneratorTreeBase_Benchmark::generateShield(TVarListHandler* xVars,
		TVarListHandler *xSupport) {
	n_shielding_queries=0;
	TShieldGeneratorTreeBase::generateShield(xVars,xSupport);
}

bool TShieldGeneratorTreeBase_Benchmark::checkCondition(int xA, int l, int yB,
		int* xMap) {
	int nIndex,xn,yn;


	for(nIndex=0;nIndex<xNeighbours->lenList->at(xA);nIndex++) {
		// extract id of neighbour
		xn=xNeighbours->varList[xA]->at(nIndex);
		yn=xMap[xn];
		n_shielding_queries++;
		if(checkConditionPlane(xA,xn,l,yB,yn)) {
			return true;
		}
	}

	return false;
}



template class TShieldGeneratorTree_SqrEuclideanPrototype<TShieldGeneratorTreeBase>;
template class TShieldGeneratorTree_SqrEuclideanPrototype<TShieldGeneratorTreeBase_Benchmark>;

template class TShieldGeneratorTree_TorusPrototype<TShieldGeneratorTreeBase>;

template class TShieldGeneratorTree_SqrEuclideanNoisePrototype<TShieldGeneratorTreeBase>;
template class TShieldGeneratorTree_SqrEuclideanNoisePrototype<TShieldGeneratorTreeBase_Benchmark>;

template class TShieldGeneratorTree_PEuclideanPrototype<TShieldGeneratorTreeBase>;
template class TShieldGeneratorTree_PEuclideanPrototype<TShieldGeneratorTreeBase_Benchmark>;

template class TShieldGeneratorTree_SpherePrototype<TShieldGeneratorTreeBase>;
template class TShieldGeneratorTree_SpherePrototype<TShieldGeneratorTreeBase_Benchmark>;


template class TShieldGeneratorTree_ReflectorPrototype<TShieldGeneratorTreeBase>;

